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Abstract: Optimal designs are required to make efficient statistical 
experiments. D-optimal designs for some models are calculated by 
using canonical moments. On the other hand, integrable systems are 
dynamical systems whose solutions can be written down concretely. 
In this paper, polynomial regression models with prior information 
are discussed. In order to calculate D-optimal designs for these mod- 
els, a useful relationship between canonical moments and discrete 
integrable systems is used. By using canonical moments and discrete 
integrable systems, an algorithm for calculating D-optimal designs 
for these models is proposed. Then some examples of applications of 
the algorithm arc introduced. 
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1 Introduction 

In statistics, design of experiments is a methodology to make efficient ex- 
periments. Optimal designs have been studied by numerous authors in the 
literature. Especially D-optimal designs have been investigated by many 
authors. One of the approaches for D-optimal designs is to use canonical 
moments. 

D-optimal designs and D s -optimal designs for polynomial regression mod- 
els were studied in [18]. In the calculation in |18| . an important point is to 
use canonical moments instead of ordinary moments, and D-optimal de- 
signs can be identified by its canonical moments. As we see in Section [21 
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the objective function is written down in term of canonical moments. Be- 
sides this, D-optimal designs for various models have been calculated in 

[2i si m ei ei ei ei nni hd mi nn eui ng eu imi- For example, d- 

optimal designs for weighted polynomial regression with a weight function 
x a (l — x)P were found explicitly in [19] by using canonical moments. How- 
ever, in most cases, an explicit form of the D-optimal design is unknown. 
We here consider polynomial regression models with some prior information. 
For example, D-optimal designs for this kind of models have been studied in 
[31 El EI]- D-optimal designs for polynomial regression models with only odd 
(or even) degree terms are calculated in [3]. Polynomial regression without 
intercept is considered in Polynomial regression models through origin 
is considered in |S]. 

On the other hand, the term integrable systems is used for nonlinear 
dynamical systems whose solutions can be written down concretely. For 
Hamiltonian systems with finite degree of freedom, their integrability is 
defined in the Liouville- Arnold theorem pp. However, even now, there is no 
mathematical definition of integrability for nonlinear systems with infinite 
degree of freedom. This is why we call an explicitly solvable nonlinear system 
as an integrable system. Discrete integrable systems are discrete analogues 
of integrable systems. That is, it is well-known that integrable systems can 
be discretized such that descretized systems are also solvable. 

Discrete integrable systems have been applied to numerical analysis. 
Typical examples are matrix eigenvalue algorithms [151 [22] , and algorithms 
to compute matrix singular values |23j . 

An algorithm for calculating D-optimal designs for polynomial regression 
through a fixed point is proposed in [17]. In |17j . the relationship between 
canonical moments and discrete integrable systems are used. In this thesis, 
we propose an algorithm for constructing D-optimal designs for polynomial 
regression with prior information, which is a generalized algorithm of [17] . 
and the considered models include polynomial regression through origin [8], 
and some weighted polynomial regression as particular cases. Moreover 
the algorithm can be applied to some optimal designs, and it allows us to 
calculate a larger class of optimal designs. That means the relationship 
between canonical moments and discrete integrable systems expands the 
class of D-optimal designs which can be calculated. 
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2 A preliminary 



This section give a brief introduction of polynomial regression models, D- 
optimal designs, and canonical moments. 

At first, we consider the following common linear regression model 

Y = 6 T f(x) + e, 
E[e] = 0, V[e] = a 2 , 

where f(x) = (fo{ x ) fi{ x ) • • • fm-i(x)) T denotes a known vector of linear 
independent functions, 6 = (6o(x) 9\{x) ■ ■ ■ # m _i(x)) T denotes an unknown 
vector of parameters, e denotes the error term. Here we assume each exper- 
iment has a stastically indendent error term e. 

Let Vi be the set of probability mesures on the borel set on /. For given 
IX G Vi, let Cfc denote the kth moment fj x k d/j,(x). Suppose the number 
of expriments be n, and experimental conditions x\,X2, ■ ■ ■ , x n should be in 
interval [0, 1]. Here we consider design x±, x%, . . . , x n is corresponding to the 
probability mesure \i G ^[o.i] such that fj,({x}) = #{k\xk = x}/n. Then D- 
optimal designs are defined as probability mesure \x € Pro i] which maximizes 
the determinant of Fisher information matrix Mf(fi) = f (x) f (x) T dfj,(x) . 
In the case of (m — l)th polynomial regression, that is, in the case where 
fk( x ) = x > the D-optimal design is defined as the optimization solution of 
the optimization problem 

maximize \a +j \™fj 
subject to fj, G ^[0,1] • 

Note that, in the optimization problem, while fi({x}) should be multiple of 
1/n for all x, we do not take a such constraint. Therefore we consider only 
the relaxation optimization problem. See |16| for details. 

In the optimization problem ([T]), the objetive function is written down 
in terms of moments. However the constraint is complicated in terms of 
moments. To simplify the constraint, sometimes the optimization problem 
is rephrased in terms of canonical moments. 

Now, we define canonical moments. Suppose we consider the set V[o,i] of 
the probability measures on [0,1]. For a given probability measure /U G ^[0,11 j 
let ct denote the maximum of the A;-th moment over the set of all measure 
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having moments cq,ci,... , cj~-i- Similarly let c k denote the corresponding 
minimum. Canonical moments are defined by 

Ck- — c, 

Pk = ^ k ~, k = l,2,...,N, (2) 

c k ~ c k 

where N is the minimum of j which satisfy Cj +l = c~ +l . If Cj > cj for an 
arbitrary positive integer j, let's set N = oo. 
Canonical moments have the property that 

U 6(0,1), fe = l,2,...,iV-l, 

We {0,1}. 

Conversely, for a given arbitrary sequence {pk\k=\ wmcn satisfies ((3|), there 
is a unique /U € ^[0,1] which has the canonical moments {pk} k =i- 

It is to be noted that canonical moments have a Hankel determinant 
expression. Let H^\h^ be the Hankel determinants defined by 

"fc ~~ l^+j+^lij^O ' ~~ \Ci+j+n L -i+3+n+l\ij = o i ^ 

(n) 

where /c = 1,2,..., re = 0, 1, 2, . . .. We here set that Hq = 1, and that 
H^ 1 = if a matrix size k is negative. The canonical moments pf. have the 
following Hankel determinant expression: 

H { k 1) H k °\ , 
P2fc " 1 " -(0)77(1) ' P2k - MDrfo) - 

We introduce useful variable Cfc by the transformation of canonical moments 
Pa- 



Co =0, Ci=Pi) 

Cfe = (1 - Pk-i)Pk, k = 2,3, ... ,7V. 
Then Cfe also have the Hankel determinant expression 

> _ -"fc -"fc-l > _ n k+l I1 k-l r, — 1 9 

- „(o) „(i) ' <> 2fc - „(1) „(0) ' ' ' 



(5) 
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From this expression, Hankel determinants H m can be expressed in a product 
of C, k or of canonical moments p k , 



m—1 

\m—k 



k=l 

m—1 \ m—1 

n (i - PuT-i-^p n (a - Ki-ifenr. 

where 2m — 2 < N. D-optimal designs for polynomial regression are cal- 
culated in [18] through the expression ©. Similar approaches for other 
D-optimal designs are described in the book [7]. 

3 An algorithm for calculating D-optimal designs 
for polynomial regression with prior information 

In this section, we consider polynomial regression with prior information. 

In Subsection 13.21 polynomial regression with prior information is de- 
fined and formulated as a linear regression. Then we proposed an algorithm 
for calculating the D-optimal design for polynomial regression with prior 
information. 

3.1 Generalized canonical moments and the nonautonomous 
discrete time Toda equation 

At first, we generalize ordinary moments. For given moments {cfc}^L , let 
ci T ^ be defined by 

(TtU{A}) _ (T) _ \ (T) (<f>) _ 
c k — c k+l AC k ' c k — c k, 

(T) 

where T denotes a multiset. And let Hf, be its Hankel determinant 
l c i+jli J=o- Then canonical moments p k and variables Ck are expressed as 

rrW rr({0,-l}) rr({0}) rr({-l}) 



P2k tt (W) rr ({-1}) ' P2k+1 rrW rr({0-l}) ' 

H k H k H k+l H k (7) 

rrW rr({0}) rr({0}) rrW [ ' 

^ rr({0}) rrW ' ^ 2k+1 rr{<t>) MM) ' 
H k H k H k+l H k 
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(T) (T) 

We define generalized canonical moments p k and variable Q by 

rr(T) rr(Tta{0 -1}) „ (Ttti{0}) rr(TW{-l}) 

(T) H k+l H k~l (T) -"fc+l -"fc 

?2fc " „-(Ta{0}) rr(TU{-l}) ' P 2fc+1 ( T ) „(TtH{0 -1}) ' 

^fc ^fe H k+l H k 

rr(T) tt(TI±>{s}) rr(TW{ S }) fr(T) 

^k - H (TU{s}) H (T) ' Wl - jj(T) H (TU{s}) ' 
k k k ~\~ 1 /c 



(8) 



It is clear that pjj^ = p/j , c\^^ = Ck ■ While generalized canonical moments 
do not satisfy like a property ([3]) , generalized canonical moments satisfy like 
a property ([5]), that is, 

AT) _ q AT) _ (T) 

Co -u, Ci -Pi , 

JT,0) n (T) \ (T) , „ , r ^ 

Q -(i-Pfc-iK > k = 2,3,...,N. 

(T) 

Additionally variables Q is a determinant solution of nonautonomous dis- 

(T) 

crete time Toda equation, namely, Q is satisfy nonautonomous discrete 
time Toda equation 



JTvlMjM) .(Tm{Ai},A a ) x _ ,(T,Ax) AT,Xi) . 
^>2fc + ( >2fc+l + A 2-C 2 fc+1 +S 2 fc+2 + A 1> 

.(T(±l{Ai},A 2 ).(Tl±){Ai},A2) _ .(T,Ai) .(T,Ai) 
t '2fc+l t '2fc+2 ~~ <>2fc+2 t '2fc+3 - 



(10) 



From (I10p . we obtain the following formula 



^(T,Ai) .(T.Ax) . _ .(T,A 2 ) ^(T,A 2 ) . 
ATM) ATM) _ ATM) ATM) 

^2k+l ^2k+2 ~ ^2k+2 ^2k+3 " 



fill 



3.2 D-optimal designs for polynomial regression with prior 
information 

We consider the (m + S — l)-th degree of polynomial regression 



m+S-l 
k=0 

E[e] = 0, V[e] = a 2 
with the S = X^'=i va l ues as prior information 



(12) 



, 0<j<l, 0<k<bj, (13) 

x=/3 j 



I. 



where g{x) = E\Y\x] = XX=*f~ 1 @kX k , bo, b\, . . . , are positive integers, 
and /3o,/3i, . . . , are arbitrary distinct real values. That means we con- 
sider the cases where the exact S values (|13p are known before experiments. 
Note that fa do not have to be in X = [0,1]. We call the model (fT2l) 
PRM m (/3,6), for short, where (3 = (fa, fa, • • • , Pi-i), & = (&o,&i, ■ ■ -M-l)- 

There are multiple ways to consider the model PRM m (/3,6) as linear 
regression models. However the D-optimal design for the model PRM m (/3, b) 
is defined uniquely, since D-optimal designs for linear regression models 
depend on a linear space spanned by bases functions only (see [7J Theorem 
5.5.1]). The D-optimal designs are formulated as the following theorem. A 
proof of the theorem is given in Appendix. 

Theorem 3.1. The D-optimal design for polynomial regression with prior 
information PRM m (/3, b) is defined as the optimal solution of the optimiza- 
tion problem 

maximize ffi?^ 

(14) 

subject to fj, € ^[0,1] • 

where the multiset T satisfies the number of fa in the T is 2b k , that is, let 
itit(x) be the multiplicity function, then 

(2b k (x = fa) 
m T (x) = < . (15) 

[ (otherwise) 

Note that the D-optimal design for PRM m (/3,6) is equivalent to the D- 
optimal design for weighted polynomial regression model 

m— 1 



Y = J2 °kx k + e, 

w{x) 



k=0 

2 '-1 

E[s]=0, V[e\ = ^— w(x) = H(x-fa)^. 

fe=0 



Here we propose an algorithm for rephrasing the optimization problem 
14]) in terms of canonical moments. In the algorithm, the two formulas (116|) . 



(|17j) and the nonautonomous discrete time Toda equation (llOj) are used. 

At first, we describe the objective function Hm ' by using the value Cj?'^ 
corresponding to generalized canonical moments. We obtain the formula 
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which is similar to © 

771—1 

^»=(4 T, )™n(«' 0, ) ro " 1 - <«» 

Then we describe the objective function H%P in terms of the value C^' ^ 
corresponding to canonical moments. Here by using the nonautonomous dis- 
crete time Toda equation (fTU|) and the formula (fTT|) . can be expressed 

in terms of ; d^' ^ > And we obtain the formula about Cq 7 ^ 

(Ta{s}) 

cf'^^r, c f = i, (it) 

c 

then Cq 7 ^ also can be expressed in terms of Ci^ , C2 > • • • by using the 
formula. 

Lastly we describe the objective function H m in terms of canonical 
moments by using the relationship ©. 

By putting it all together, the proposed algorithm for calculating the 
D-optimal design for polynomial regression with prior information (|12p is 
described as follows. 



The algorithm for calculating D-optimal designs for 

PRM m (/?,&) 



Step 1. By using the formula (|T6l) . describe the objective function 
Hm in terms of (j?'^ and . 

Step 2. By using the nonautonomous discrete time Toda equation 
(fTUj) and the formulas (fTT|) . (|T7|) . describe the objective func- 
tion fl^ in terms of Cfc^'°^ = Ck- 

Step 3. By using the relationship ([5]), describe the objective function 

(T) 

H m in terms of canonical moments 
Step 4. Find canonical moments which maximize the objective func- 

(T) 

tion H m . 



4 Application of the algorithm for calculating some 
D-optimal designs 

In this section, we introduce two examples of application of our algorithm. 
In Subsection 14. 1\ robust D-optimal designs for approximate polynomial 
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regression with prior information is considered. This model is generalized 
model considered by |10j . In Subsection 14.21 maximin optimal designs for 
estimating a function of parameters on weighted polynomial regression. We 
consider the generalized case of the case considered by [6]. 

4.1 Robust D-optimal designs for approximate polynomial 
regression with prior information 

In this subsection, let the design space be X = [—1, 1] instead of [0, 1]. There 
is one to one correspondence between [i G ^[0,1] and a symmetric measure 
£ G ^[-1,1] such that 

The approximate polynomial regression model is described as 

m—l 

Y=Y, QkX k + x m ^(x) + e, 

k=0 

E[e] = 0, V[e) = a 2 , 

where ip denotes an unknown function. Let the best linear unbiased es- 
timator be 0{ip) when we estimate as ip(x) is considered as 0. Then we 
obtain 

B[^)-fl(0)]=5(f)(0^) 

where £ is a design, n is the number of observations, and 
r(V>) = J (1, x, . . . , x m - 1 ) T x m ^(x)d^(x), 

For given a continuous function r](x) and a positive number d, maximin 
optimal design is defined as 

maximize Hff (£ ) 

subject to £ G V hlA] , sup r(^) T (B$ )" 1 r(V') < d. 
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In [10], the case where 

r)(x) = \x\ a , a e Z> (18) 

is considered, and it is shown that the constraint sup^i^i r(^) T (5m^)" 1 r(?/)) < 
d is described in terms of canonical moments as 

L(m+a)/2j m+a-2i 

sup r^fiB^r'rW = £ S i>m+a ^f [] Cj(m), (19) 
I^M i=K2j+i j=i 

where Sij(/j,) is defined recursively as 

S id (n) = 0, 0<j<i, 
Sij(fi) = 1, i = 0, i > 0, 

Sijifi) = + Cj-i+i^Si-ijG"), < i < j. 

Now we consider the approximate polynomial regression model with prior 
information, namely, 



25+m-l 

Y= °k xk + x 2S+m ^{x) + e, 

k=0 

E[e] = 0, V[e] = a 2 , 



(20) 



with the 2S values as prior information 



dx^ 9[X) 
d k 

dx~ k9{x) 



x=/3 j 

, 0<j<l, 0<k<bj, 



(21) 



X=-Pj 

where g(x) = E[Y\x] = Y^T^a" 1 ®kX k , and /3o, /Si, . . . , fii-i are arbitrary 
distinct real values. Note that prior information (|2ip must be symmetric 
with respect to the origin. 

We can show that the optimization problem corresponding to the model 
(|20D is the following by similar calculation to |10j 



maximize 

subject to £ € sup r( T ')(^) T ( J B^'))~ 1 r( T ')(V') < d, ^ 
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where T' denotes the multiset satisfying 
m,T'(x) - 



2b k {x = /3 k ) 
2b k (x = -/3 k ) 
(otherwise) 



and 

>(T')fc\ - (. 

-i+j H,j=0 



r (T ' } (V) = J\l,x, . . ,x m4 ) T x"VW \J[(x - Pi)(x + foj d£(x). 

The constraint of (|22p corresponding to (|19p is described in terms of 
generalized canonical moments as 

L(m+a)/2j m+a -2i 

(cTco)- 1 £ sgU*(M) 2 n cf°V), 

i=La/2j+l j=l 

where T is the same as (flBT) . and 



5g ) ( / u) = 0, 0<j<*, 
SgV) = l, i = 0, i>0, 

CV) = C-i(^) + C^2i(M)5Sj(M), < i < j. 

Hence we can obtain the expression of the optimization problem correspond- 
ing (|20p in terms of canonical moments by using our algorithm. 

4.2 Maximin optimal designs for estimating a function of 
parameters on weighted polynomial regression 

Consider the weighted polynomial regression 

m—1 



Y = J2 9 h* k + e, 

k=0 

l-l 

E[e] = 0, V[e] = a 2 w(x), w(x) = J[(x- p k ? bk • 



k=0 
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In this subsection, we consider an optimal design for estimating ff m -l(^m-l)+ 
9m-2(8m-2) + ■ ■ ■ + go(#o)> where g k is a polynomial. In this case, the inverse 
of the asymptotic variance of the estimator is 



m-l , , \ 2 

fc=0 ^ fc ' 



where ip k (n) = Hj,^/ H k , and T is the same as (fT5]h Since the variance 
of estimator depends on the unknown parameters 6^, we consider maximin 
optimal design defined as the optimization solution of the optimization prob- 
lem 

maximize min7(u, 9) 

eee 

subject to fi S ^[0,1] 

where is a given parameter space. Suppose the parameter space © = 
{0\s k <6 k < t k }. 

Prom [6l Theorem 3.1], the optimal solution of the optimization problem 

maximize / 7(/i, #) p d7r(#) 

Je (23) 

subject to /i € "P[o,i]- 
converges weakly to the maximin optimal design as p — >• — oo, where 

m— 1 



d7r(0) = d6 n h k (9 k ), 
hk(0k) = 



fc=0 

2 



1 (otherwise) 

Then we can calculate the integral in the objective function of the (|23p . 
And we can express in terms of canonical moments by our algo- 

rithm. Therefore, after expressing the optimization problem (|23p in terms 
of canonical moments, we can obtain an approximate maximin optimal de- 
sign by solving (j23[) numerically for small p. 
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A The proof of Theorem 13.11 

It can be turn out that the optimization problem (|14p corresponding to 
PRM m (/3, b) corresponds to the vector of basis functions 

/(*)= fn^-ft) 6j j (i,^-.-^ m - i ) T - (24) 

Let gk(x) = (1, x, . . . , x fc_1 ) be the vector of basis functions for polyno- 
mial regression, then the vector (|24|) of basis functions corresponding to 
PRM m ((/3 ,/3i, • • • , A-i), Q>o,bi, ■ ■ ■ ,h-i)) is expressed as 

l-l 

f(x) = H(x-l3 j pg m (x). (25) 

5=0 

To prove Theorem 13. 1\ it suffices to show that PRM m ((/?o, • • • , A-i)> 
(6q, &i, . . . , &z-i)) corresponds to the vector ([25]) of basis functions. Let us 
prove it by the principle of induction. By the symmetricity of (3j and bj, 
we can assume that PRM m +i((/3o, Pi, • • • , Pi-i), (bo, b± — l,..., 6/_i)) corre- 
sponds to the vector of basis functions 

/(*) = (x- Pof ' 1 m (x - Pj) h A g m+ i(x). 

Let M(x) = Y^j=i( x ~Pj) bj i then the linear regression model PRM m+1 ((/3 , /3 l5 . . . 
(6o — 1, 6i, . . . , h-i)) is described as 

m 

Y = (x — po) b °~ l M(x) Qkx k + e. 

k=0 

Thus the linear regression model PRM m ((/3o, Pi, ... , Pi-i), (bo, b±, . . . , fy-i)) 
is described as 

m 

Y = (x- Po) b °~ 1 M(x) k x k + e. (26) 

fc=0 

with one given value as prior information 
d b °~ 



dx b o 



-A{x-Po) h °- l M{x)Y,hx k 



k=0 



(27) 



X=fio 
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Since 



dx k 

the value (1271) becomes 



(x-(3 ) b °- 1 



0, k = 0,1,..., b 



x=p 



k=0 



by the general Leibniz rule. Hence we obtain the value 



(28) 



k=0 



from prior information ()27|) . 

Prom (|28|) . substitute 9q = a — SfcLi ^fc/^o hito the model (f26|) . we obtain 

(mm \ 
J^^-^^ + al+e. (29) 
fe=i fc=i / 

When we obtain a response by the obserbation at the experimental con- 
dition Xfc, we can calculate the value yk — {%k — Po) b °~ 1 M(xk)a easily. Thus 
we can ignore the term (x — (3o) bo ~ l M (x)a from the model (I29f) . then we 
obtain the model 



Y = (x — Po) b °- l M{x) J2 Ux k ~ Po) + e. 



(30) 



k=l 



Here the vector of basis functions corresponding to the model fj30f) is 

f(x) = (x- f3 ) b °- l M(x) 



x 2 - ft 



„m ora 



(31) 



Let the non-singular matrix A be 





( 1 





•• 


• o\ 




-/3o 


1 


•• 


• 


A = 





-/3o 


1 •• 


• 




I 





•• 


• V 
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then, by multiplying A to the vector (|3ip of basis functions from left, 

/ x-Pq \ 

x 2 - ffi - p (x - /3 ) 



Af(x) = {x- fo) b °- l M{x) 



(x - /3 ) fe °" 1 M(x) 



y x m_pm_^ x m-l_ p m-l ) j 

( x-p \ 
x(x - Pq) 



\x m - 1 (x-p )J 
(x-p ) b °M(x)g m (x). 



Thus we obtain ([25lh 
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